class: inverse,left, middle background-image: url(data:image/png;base64,#background.png) background-size: cover <img src="data:image/png;base64,#LOGO_DIPLOMADO.png" width="500px"/> ##MODULO 1: Manipulación y análisis de datos geoespaciales ### Analisis Espacial Gabriel Castro. B | Javiera Aguayo .T <br> Correos: gabriel.castro.b@pucv.cl | javiera.aguayo@pucv.cl .large[<b><a href="mailto:xxxxx@pucv.cl">LabGRS</a> | Agosto, 2022</b>] <br> --- class: center,middle background-image: url(data:image/png;base64,#labgrs_logo.png) background-size: 35% --- ##Contenidos **Analisis Espacial** **Tipos de Operaciones:** **a) Operaciones de Consulta** **b) Operaciones de Transformación** **c) Operaciones de Alteración** **d) Otras Operaciones Complementarias** --- ## Análisis Espacial .pull-left[ Es un análisis geográfico que tiene como objetivo explicar patrones o relaciones espaciales en base a términos u operaciones matemáticas, y esto vincula diectamente a la aplicación de operaciones en objetos tipo raster, ya que cada celda nos indica información espacial relevante en términos de localización. ] <right><img src="data:image/png;base64,#dibujo_raster.png" width="500px"/></right> --- .pull-left[ ## Tipos de Operaciones Espaciales Al igual que con los objetos tipo vector, con los objetos tipo raster también se pueden llevar a cabo una serie de operaciones basicas de las mismas tres tipologías anteriores: **Consultas**:Se utiliza para revisar las características del objeto raster (Sistema de proyección, extensión, n°de bandas, etc.) **Transformación** :Operaciones que modifican en la geo-referencia o el tipo de estructura del objeto raster. **Alteraciones**: Son operaciones que modifican tanto la información temática como también la geométrica del objeto raster. ] -- .pull-right[ <right><img src="data:image/png;base64,#tipos_de_procesos.png" width="500px"/></right> ] --- ## Diagrama de operaciones <center><img src="data:image/png;base64,#Diagrama_operaciones_analisis_espacial.png" width="700px"/></center> --- Librerías utilizadas en esta sesión: ```r library(terra) library(magrittr) ``` Objeto a utilizar en esta sesión ```r DEM_01<-rast('s35w072.hgt')# Modelo elevación digital. DEM_01 ``` ``` ## class : SpatRaster ## dimensions : 3601, 3601, 1 (nrow, ncol, nlyr) ## resolution : 0.0002777778, 0.0002777778 (x, y) ## extent : -72.00014, -70.99986, -35.00014, -33.99986 (xmin, xmax, ymin, ymax) ## coord. ref. : lon/lat WGS 84 (EPSG:4326) ## source : s35w072.hgt ## name : s35w072 ``` --- ### Operaciones de Consulta: **Clase (class):** - Indica la tipología del objeto que estamos abriendo ```r class(DEM_01) ``` ``` ## [1] "SpatRaster" ## attr(,"package") ## [1] "terra" ``` En este caso la clase del objeto es un "SpatRaster", que es un tipo de objeto raster correspondiente al paquete Terra, el cual se puede componer de una o multiples capas. **Sistema de Proyección (crs):** Muestra el sistema de coordenadas o tipo de proyección del objeto raster. ```r crs(DEM_01, proj=T, describe=T, parse=F) ``` ``` ## name authority code area extent proj ## 1 WGS 84 EPSG 4326 <NA> NA, NA, NA, NA +proj=longlat +datum=WGS84 +no_defs ``` --- ### Operaciones de Consulta **Histogramas (hist):** #### Los histogramas nos permiten... .pull-left[ - **Revisar la cantidad de pixeles (frecuencia) de nuestro archivo Raster para un valor especifico** - **Revisar la distribución de los pixeles en funcion de los valores que almacenan** - **Para el caso de un DEM se observa la frecuencia de pixeles segun la elevación** ] .pull-right[ ```r Histograma_dem <- terra::hist(DEM_01) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-5-1.png" width="100%" /> ] --- ####Podemos mejorar la representacion de la información a traves de los elementos de la funcion hist(): ```r histograma_mejorado <- terra::hist(DEM_01,maxcell=ncell(DEM_01), plot=TRUE, main = "Histograma DEM región Maule", col="blue", xlab = "Elevación", ylab = "Frecuencia (N°pixeles)", cex.lab=0.8, cex.axis=0.8, cex.main=0.8, cex.sub=0.8) # Ajustar parametros según la visulización deseada ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-6-1.png" width="100%" /> --- ```r class(histograma_mejorado) ``` ``` ## [1] "histogram" ``` - ***Objeto creado es de tipo "histograma"*** ```r print(histograma_mejorado$counts) # frecuencia de pixeles por rango de elevación ``` ``` ## [1] 139489 938895 4498624 3392956 1863107 827987 545829 329036 156096 ## [10] 95005 67850 46298 31027 14467 8226 4858 3457 2384 ## [19] 1245 229 136 ``` ```r print(histograma_mejorado$mids) # volores intermedios del rango de elevación ``` ``` ## [1] -50 50 150 250 350 450 550 650 750 850 950 1050 1150 1250 1350 ## [16] 1450 1550 1650 1750 1850 1950 ``` -- **- Los parametros de nuestro histograma como los intervalos de elevación y el rango altitudinal pueden ser modificados por el usuario** **- Esto ayuda a consultar información especifica de nuestro raster** --- ### ¿Que pasa si tengo muchos archivos rasters? ¿Es posible graficar todo de una vez? #### Ej: Se tiene un total de 3 DEM´s a los cuales quiero consultar por histograma: ```r DEM_01 <- rast("s35w072.hgt") DEM_02 <- rast("s37w073.hgt") DEM_03 <- rast("s36w072.hgt") lista_DEM <- list(DEM_01,DEM_02,DEM_03) print(lista_DEM[1])# ver el primer elemento de la lista ``` ``` ## [[1]] ## class : SpatRaster ## dimensions : 3601, 3601, 1 (nrow, ncol, nlyr) ## resolution : 0.0002777778, 0.0002777778 (x, y) ## extent : -72.00014, -70.99986, -35.00014, -33.99986 (xmin, xmax, ymin, ymax) ## coord. ref. : lon/lat WGS 84 (EPSG:4326) ## source : s35w072.hgt ## name : s35w072 ``` --- #### Aplicar ciclo for a nuestra lista para obtener los 3 histogramas. Esto ayuda a la optimizacion de tiempo y codigo. ```r a<-list() op <- par(mfrow=c(1, 3)) for (i in 1:length(lista_DEM)) { a[[i]] <- hist(lista_DEM[[i]], maxcell=10e10, plot=TRUE, ain = "Histograma DEM región Maule", col="blue", xlab = "Elevación", ylab = "Frecuencia (N°pixeles)") } ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-10-1.png" width="100%" /> ```r par(op) ``` --- ### Operaciones de Consulta **Resumen (summary):** Muestra un resumen de las principales estadísticas del contenido del objeto, en este caso de tipo ráster. ```r summary(DEM_01, size=1e20) # argumento size nos ayuda a tomar toda la informacion del ráster. ``` ``` ## s35w072 ## Min. : -20.0 ## 1st Qu.: 148.0 ## Median : 222.0 ## Mean : 269.4 ## 3rd Qu.: 330.0 ## Max. :1969.0 ``` - Elevación máxima: 1969 m. - Elevación mínima: -20 m. - Elevación promedio: 269 m. Recordar que es DEM es de la zona costera de Constitución --- ## Operaciones de Transformación **Reproyectar (project):** Cambia el sistema de coordenadas de un objeto raster. ```r DEM_reproyectado <- project(DEM_01, "epsg:32719") #cambio de proyeccion A WGS84 / UTM HUSO 19S DEM_reproyectado ``` ``` ## class : SpatRaster ## dimensions : 4002, 3343, 1 (nrow, ncol, nlyr) ## resolution : 28.29421, 28.29421 (x, y) ## extent : 222895.4, 317483, 6122823, 6236057 (xmin, xmax, ymin, ymax) ## coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) ## source : memory ## name : s35w072 ## min value : -18.5506 ## max value : 1967.236 ``` --- ### Operaciones de Transformación **Reproyectar (project):** .pull-left[ ```r plot(DEM_01, main = "DEM coords geograficas") ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-13-1.png" width="100%" /> ] .pull-right[ ```r plot(DEM_reproyectado, main = "DEM coords UTM") ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-14-1.png" width="100%" /> ] --- ### Operaciones de Transformación **Guardar raster (writeRaster):** Guarda un nuevo objeto tipo raster ```r writeRaster(DEM_reproyectado, filename = "DEM_UTM.tif", #Nombre para guardar el objeto en disco datatype = "FLTS4", #Formato del tipo de dato con el que se guardará el objeto overwrite = T) ``` --- Si queremos revisar el tipo de data type que tiene un objeto, podemos utilizar las siguientes líneas: ```r DEM_reproyectado %>% raster::raster() %>% raster::dataType() ``` ``` ## [1] "FLT4S" ``` <right><img src="data:image/png;base64,#data_type.png" width="500px"/></right> https://www.rdocumentation.org/ --- ### Operaciones de Transformación **Vectorización (as.polygons):** Convierte un objeto tipo raster en un objeto tipo vector. En este caso en específico en un SpatVector de tipo polígono. ```r Land_Cover<-rast('Land_cover.tif',) plot(Land_Cover, col = terrain.colors(7), #Paleta de colores type="classes") #Tipo de leyenda ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-17-1.png" width="100%" /> --- ```r poli_LC<-as.polygons(Land_Cover) #Vectorización de un raster de usos de suelo plot(poli_LC,"LC_2014_Zhao", col=rainbow(7), #Paleta de colores border="black", #Borde de lineas lwd=1, # Grosor de línea type="classes", # Tipo de leyenda legend=T ) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-18-1.png" width="100%" /> La vectorización puede hacerse con diferentes funciones, dependiento de la tipología del objeto (polígono, línea o punto): ### -póligono -> **as.polygons()** ### -línea -> **as.lines()** ### -puntos -> **as.points()** ### -contorno -> **as.contour()** --- ### Operaciones de Alteración **Unión (merge):** Se utiliza cuando se quiere formar un nuevo objeto raster con una mayor extenxión espacial, a partir de varios objetos raster. ```r # Crear una lista de los objetos raster DEM_list<-list.files(path="~/LABGRS/Diplomado/MDO01_Clase_de_Analisis_Espacial", #Directorio de origen pattern=glob2rx("s*.hgt"),# full.names=T, recursive=T) lista<-list() # Lista vacía for(i in 1:length(DEM_list)){ lista[[i]]<-DEM_list[i] %>% rast } lista[[1]] ``` ``` ## class : SpatRaster ## dimensions : 3601, 3601, 1 (nrow, ncol, nlyr) ## resolution : 0.0002777778, 0.0002777778 (x, y) ## extent : -72.00014, -70.99986, -35.00014, -33.99986 (xmin, xmax, ymin, ymax) ## coord. ref. : lon/lat WGS 84 (EPSG:4326) ## source : s35w072.hgt ## name : s35w072 ``` --- Una vez creada la lista de objetos raster debemos utilizar la función **sprc()** el cual nosrealizar un collección de raster sin restricción de que los rasters sean similares entre sí. Y finalmentela función **merge()** ```r DEM_Merge<-lista %>% #llamar lista de objetos raster terra::sprc() %>% # Crear un objeto tipo SpatRasterCollection merge() #Unir todos los objetos del SpatRasterCollection plot(DEM_Merge) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-20-1.png" width="100%" /> --- ### Operaciones de Alteración **Recorte (crop):** Se utiliza para realizar un recorte de la extención de un objeto raster en base a un objeto tipo vector. Para esto se debe abrir un objeto tipo vector con la función **vect()**, y posteriormente podemos hacer un reproyección entre los objetos para que estos tengan la misma proyección. ```r #Recorte prov_cauquenes<-vect("Provincia_Cauquenes.shp") DEM_UTM<-project(DEM_Merge, "epsg:32719", method='near', res=30 ) DEM_cr<-crop(DEM_UTM, prov_cauquenes) ``` --- ```r plot(DEM_cr) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-22-1.png" width="100%" /> --- ### Operaciones de Alteración **Máscara (mask):** Se aplica una máscara a los valores que se encuentren dentro de un área específica, que puede ser indicada por un objeto tipo vector. Aquellas celdas con valores que no se encuentren dentro del área, se les asignará un nuevo valor (NA de forma predeterminada). ```r #Máscara DEM_cauquenes<-mask(DEM_cr, prov_cauquenes, inverse=F, #La máscara se aplica aplica a los valores fuera los límites del área de interés. touches=T) #La máscara no descarta aquellos valores sobre el límite del polígono ``` .pull-left[ ```r plot(DEM_cauquenes) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-24-1.png" width="100%" /> ] <rigth><img src="data:image/png;base64,#touches_mask.png" width="400px"/></rigth> --- ### Operaciones de Alteración **Reclasificación (Classify):** Corresponde a la asignación de nuevo valores a un dato o conjuntos de datos. Utilizaremos la función **classify()**, donde se puede realizar la clasificación en base a un obejeto tipo matríz de tres columnas. ```r m<-c(-Inf,250,1,250,500,2,500,750,3,750,Inf) m ``` ``` ## [1] -Inf 250 1 250 500 2 500 750 3 750 Inf ``` ```r mt_DEM<-matrix(m, nrow=4, ncol=3, byrow = T ) mt_DEM ``` ``` ## [,1] [,2] [,3] ## [1,] -Inf 250 1 ## [2,] 250 500 2 ## [3,] 500 750 3 ## [4,] 750 Inf -Inf ``` --- ```r class_DEM<-classify(DEM_cauquenes, mt_DEM) class_DEM ``` ``` ## class : SpatRaster ## dimensions : 2647, 2413, 1 (nrow, ncol, nlyr) ## resolution : 30, 30 (x, y) ## extent : 158622, 231012, 5982103, 6061513 (xmin, xmax, ymin, ymax) ## coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) ## source : memory ## name : s35w072 ## min value : - ? ## max value : 3 ``` ```r plot(class_DEM) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-27-1.png" width="100%" /> --- ### Operaciones de Alteración **Resampling (resample):** Se realiza un resampling o remuestreo de datos cuando dos raster que se quieren homologar o relacionar, tienen caracterísitcas geometricas distitas, tales como: sistema de coordenadas, extensión y/o resolución espacial. Existen diferentes tipos de métodos para la estimación de los valores de las celdas, tales como: vecino mas cercano, interpolación bilineal, interpolación cúbica, etc. ```r #Resampling r_vacio<-rast(extent=ext(DEM_cauquenes), resolution=100, crs= "epsg:32719") r_vacio ``` ``` ## class : SpatRaster ## dimensions : 794, 724, 1 (nrow, ncol, nlyr) ## resolution : 100, 100 (x, y) ## extent : 158622, 231022, 5982103, 6061503 (xmin, xmax, ymin, ymax) ## coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) ``` --- ```r DEM_100<-resample(DEM_cauquenes, r_vacio, method="near") DEM_100 ``` ``` ## class : SpatRaster ## dimensions : 794, 724, 1 (nrow, ncol, nlyr) ## resolution : 100, 100 (x, y) ## extent : 158622, 231022, 5982103, 6061503 (xmin, xmax, ymin, ymax) ## coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) ## source : memory ## name : s35w072 ## min value : -3 ## max value : 802 ``` ```r plot(DEM_100) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-30-1.png" width="100%" /> --- ### Otras operaciones complementarias **1. Operaciones para extración de valores:** #### 1.1. Extracción de valores en base a un objeto tipo vector **Pregunta:** ¿Cuáles son los valores de píxeles, con la altura máxima, en las tres comunas de la Provincia de Cauquenes? -- ```r altura_max <- terra::extract(x=DEM_cauquenes, y=prov_cauquenes, fun=max, na.rm=T, touches=T) #weights = para valor ponderado as.data.frame(altura_max) ``` ``` ## ID s35w072 ## 1 1 805 ## 2 2 731 ## 3 3 801 ``` --- #### 1.2. Extracción de valores de varios pixeles ```r valores_pix <- values(DEM_cauquenes) %>% na.omit %>% # Se omiten aquellos pixeles NA as.numeric() # extraemos todos los valores de todos los pixeles head(valores_pix) ``` ``` ## [1] 160 164 160 163 175 177 ``` -- #### 1.3.Extracción de un valor de un pixel en expecifico segun su xy ```r plot(DEM_cauquenes) xy <- click() #hacer clic en el plot para extraer coords p <- cellFromXY(DEM_cauquenes, xy) #obtenemos la posicion de la celda en el DEM_cauquenes segun nuestro xy pix_val <- as.numeric(DEM_cauquenes[p]) pix_val ``` --- **2. Cálculo de área:** **Pregunta:** ¿Cuánta superficie del DEM, tiene alturas superiores o iguales a 600 msnm? -- ```r rDEM <- DEM_cauquenes # Darle un nuevo nombre al objeto DEM rDEM ``` ``` ## class : SpatRaster ## dimensions : 2647, 2413, 1 (nrow, ncol, nlyr) ## resolution : 30, 30 (x, y) ## extent : 158622, 231012, 5982103, 6061513 (xmin, xmax, ymin, ymax) ## coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) ## source : memory ## name : s35w072 ## min value : -4 ## max value : 805 ``` --- **Pregunta:** ¿Cuánta superficie del DEM, tiene alturas superiores o iguales a 600 msnm? ```r rDEM[rDEM <= 600] <- NA # Valores <= 600 msnm, sean clasificados como NA rDEM[!is.na(rDEM)] <- 1 # Valores que no son NA, sean clasificados con valor 1 rDEM ``` ``` ## class : SpatRaster ## dimensions : 2647, 2413, 1 (nrow, ncol, nlyr) ## resolution : 30, 30 (x, y) ## extent : 158622, 231012, 5982103, 6061513 (xmin, xmax, ymin, ymax) ## coord. ref. : WGS 84 / UTM zone 19S (EPSG:32719) ## source : memory ## name : s35w072 ## min value : 1 ## max value : 1 ``` --- ```r plot(rDEM,col="red") ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-36-1.png" width="100%" /> ```r total_pixeles <- values(rDEM) %>% na.omit() %>% sum() superficie<-total_pixeles*30*30/1000000 #Total de pixeles con valor 1, por el tamaño del pixel (30*30), partido por 1000000, para que quede en km2 superficie ``` ``` ## [1] 56.2743 ``` --- **3. Transformación de un Objeto SpatVector a un sf:** **Pregunta:** ¿Cómo tranformo un objeto SpatVector a un sf? -- ```r cauquenes_sf <- sf::st_as_sf(prov_cauquenes) # transformamos la clase de objeto de SpatVector a sf class(cauquenes_sf) ``` ``` ## [1] "sf" "data.frame" ``` ```r as.data.frame(cauquenes_sf) ``` ``` ## CUT_REG CUT_PROV CUT_COM REGION PROVINCIA COMUNA SUPERFICIE ## 1 07 072 07201 Maule Cauquenes Cauquenes 2126.16 ## 2 07 072 07202 Maule Cauquenes Chanco 528.01 ## 3 07 072 07203 Maule Cauquenes Pelluhue 369.56 ## geometry ## 1 POLYGON ((200501 6044948, 2... ## 2 POLYGON ((188523.5 6061237,... ## 3 POLYGON ((178639.4 6035069,... ``` --- **4. Agregar una nueva entidad espacial a un objeto tipo vector:** **Pregunta:** ¿Cuál es la altura máxima de cada una de las comunas de la provincia de Cauquenes? -- ```r cauquenes_sf$ALTURA_MAX <- altura_max[,2]# agregamos a nuestro vector la info nueva cauquenes_sf ``` ``` ## Simple feature collection with 3 features and 8 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: 158631.7 ymin: 5982112 xmax: 231004.9 ymax: 6061506 ## Projected CRS: WGS 84 / UTM zone 19S ## CUT_REG CUT_PROV CUT_COM REGION PROVINCIA COMUNA SUPERFICIE ## 1 07 072 07201 Maule Cauquenes Cauquenes 2126.16 ## 2 07 072 07202 Maule Cauquenes Chanco 528.01 ## 3 07 072 07203 Maule Cauquenes Pelluhue 369.56 ## geometry ALTURA_MAX ## 1 POLYGON ((200501 6044948, 2... 805 ## 2 POLYGON ((188523.5 6061237,... 731 ## 3 POLYGON ((178639.4 6035069,... 801 ``` --- **Pregunta:** ¿Cuál es la altura máxima de cada una de las comunas de la provincia de Cauquenes? ```r plot(cauquenes_sf['ALTURA_MAX']) ``` <img src="data:image/png;base64,#Mod_01_Clase_analisis_espacial_files/figure-html/unnamed-chunk-39-1.png" width="100%" /> --- ### Bibliografía -2022. Jerry Davis. Introduction to Environmental Data Science. Capitulo 8.Raster Spatial Analysis. https://bookdown.org/igisc/EnvDataSci/raster-spatial-analysis.html -2022. Robert J. Hijmans. Terra: Spatial Data Analysis. https://rspatial.org/terra/ -2022. Robert J. Hijmans. terra.pdf. https://cran.r-project.org/web/packages/terra/terra.pdf -2012. Essentials of Geographic Information Systems. Chapter 8.Geospatial Analysis II: Raster Data. https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/s12-geospatial-analysis-ii-raster-.html --- class: inverse middle